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Abstract 

Numerical accuracy for the boundary and domain methods of the material derivative 
approach to shape design sensitivities is investigated through the use of mesh refinement. 
The results show that the domain method is generally more accurate than the boundary 
method, using the finite element technique. It is also shown that the domain method is 
equivalent, under certain assumptions, to the implicit differentiation approach not only 
theoretically but also numerically. 


Introduction 

Haug and Choi et al. 1-4 developed a unified theory of structural design sensitivity 
analysis for linear elastic structures, using a variational formulation of the structural state 
equations. This theory allows one to take the total derivative, or material derivative, of 
the variational state equation and to use an adjoint variable technique for design sensitiv- 
ity analysis. The main attraction of this approach is that it allows one to compute the 
derivatives of structural performances analytically. No discretization approximations are 
involved during the derivation, and a step size need not be specified in the calculation. 
However, the formulation requires evaluating accurate stress quantities on the boundaries 
which are often difficult to obtain. 

Accuracy of the shape design sensitivity theory was studied in Ref. 5 through the 
equilibrium condition for different types of finite elements. However, a systematic study 
through the refinement of the finite element mesh was still not found in the literature. 

To improve the accuracy of shape design sensitivities. Choi et al. 6 proposed a new 
domain method which transforms the boundary integrals into domain integrals and there- 
fore is less influenced by the the inaccurate boundary stress evaluation. This method takes 
advantage of the averaging nature of the finite element method, and is found to be more 
accurate than the boundary approach which evaluates the derivatives using boundary in- 
formation only 1-3 . However, a velocity field for the physical domain needs to be defined. 
The necessity of defining the domain velocity may indicate that this method is closely 
related to the implicit differentiation approach which also requires knowledge about the 
domain change for differentiation of the elemental stiffness matrix ". 

In this paper, the accuracy of the design sensitivity is studied through finite element 
mesh refinement for a cantilever thin plate. Results of the domain and boundary methods 
for the material derivative approach and the implicit differentiation approach are shown 
and compared. 

In a previous paper 8 , the boundary integral formulation was shown to be equivalent 
to the implicit differentiation approach, theoretically. In this report, the domain method 
is shown to be equivalent to the implicit differentiation method, both in theoretical and 
numerical aspects. 


Shape Design Sensitivity Analysis 

Two approaches for shape design sensitivity analysis are found in the literature. One is 
the well known implicit differentiation approach and the other is the variational or material 
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derivative approach. Detailed information for these two approaches is available in Refs. 4, 
7, and 8. Only a brief background is provided in the following. 

For the implicit differentiation approach, the displacement sensitivity is obtained by 
differentiating the discretized structural system of equations 

Kz — F (1) 


By assuming that the force vector F is independent of design, this leads to 

(2) 

where K is the global stiffness matrix, z is the displacement vector, and bj is the design 
variable. 

The variational design sensitivity theory uses the material derivative concept of con- 
tinuum mechanics and an adjoint variable method to obtain computable expressions for 
the effect of shape variation on the functionals arising in the shape design problem. The 
variation of the displacement functional xp with respect to shape change is derived by differ- 
entiating the variational equilibrium equation and employing the adjoint variable method, 
to obtain 1-4 

SxP = J o' 3 {z)t ij {\)V k n k dT (3) 

where xp is defined by 



in which x is the point of interest, 6 the Dirac-measure at zero, O the physical domain, a tJ 
and e tJ the stress and strain tensors, respectively. V the design perturbation and can be 
thought of as velocity, and n k the unit normal vector of moving boundary T. The vectors 
2 and A are the displacement vectors for state and adjoint equations, respectively, which 
can be expressed as follows: 


r\ i -ft- r\1 2 


db. 


db t 


= I T t z t dT 

(5) 

J r 2 

/ 6{x-x)JdU 

(6) 


where Tj is a traction vector, T 2 a loaded boundary, and - indicates the virtual displace- 
ments that satisfy the kinematically admissible displacement field. The Einstein summa- 
tion convention for a repeated index is used throughout this paper. 

To obtain Eq. 3, the traction vector T,, the kinematically constrained boundary, and 
the loaded boundary are assumed to be fixed during the design process, i.e., they are 
independent of design. Note that the variation of the displacement functional of Eq. 3 is 
only affected by the normal movement of the boundary of the physical domain. 

Physically, the adjoint equation of Eq. 6 is interpreted by applying a unit load at the 
point x , where the displacement sensitivity is of interest. In Eq. 3, one sees that only 
the boundary stress information is needed for evaluating the variation of the displacement 
functional. Unfortunately, finite element analysis usually does not provide high quality 
stress results, especially on the boundary. It has been shown that better finite element 
results give better design sensitivity estimates, by examining the equilibrium equations for 
different finite elements 5 . 
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Domain Method 


The basic idea for the domain integral method is to take advantage of the averaging 
nature of the finite element method, instead of evaluating the less accurate stresses on the 
boundary. Since the finite element method is well known to provide better solutions inside 
the finite element, the domain method has the advantage of predicting better sensitivities. 

Applying the same procedure as in obtaining Eq. 3 with the domain method, the first 
variation of the displacement constraint functional of Eq. 4, is obtained as 4,6 

**=■ J n + <« ( 7 ) 

One should note that Eq. 7 is more general than Eq. 3, since only the loaded boundary 
is assumed to be fixed, while both loaded and kinematically constrained boundaries are 
assumed to be unchanged in Eq. 3. The kinematically constrained boundary and interface 
boundary terms appear when the divergence theorem is used to transform the domain 
integral to the boundary integral. It was shown in Ref. 6 that for an interface or built-up 
structure problem this method simplifies the formulation and avoids specifying tedious 
interface conditions and provides increased accuracy for shape design sensitivities. 

To have a better understanding of Eq. 7, each term will be examined individually. 
First, since the stress tensor a* 3 is symmetric, the first term of Eq. 7 is divided into two 
parts and then integrated by parts to obtain 

/ n c i3 (X)z uk v ktJ dn = | n [z t , k v k ^ + z hk v kti \ dn 

= f T [z t ,kV k n 3 + z hk V k n t ] dT (8) 

- X )\ hh + dU 

By assuming that only the free boundary is varied, the first term of Eq. 8 disappears, 
since the traction vector is zero, i.e., a l] {\)n ] = o t3 (\)n t = 0. The second term of Eq. 8 
then can be further modified to 

L + *>>■] Vt in = /„ + *>.*] dn 

= / n S’W eS(*)V* dn 

where the velocity V*. can be parametrized as (dx k /db t )6b t . in which x k is the position 
vector and b t is the design variable. Since all the integrals are linear in design, one can 
eliminate Sb t or choose the value as a unit number. By interpreting the adjoint variable 
A as the inverse of the reduced stiffness matrix K if all the displacement sensitivities are 
desired, Eq. 9 is discretized, using the finite element formulation, as 8 
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where the subscript i with a prime superscript indicates the derivative with respect to 
the i th design variable. K, D, and B represent the stiffness, elasticity, and strain recovery 
matrices, respectively. The stress-displacement and strain-displacement relationships are 
employed in obtaining Eq. 10, which are defined in the following 

<r ,J (A) = DBX 
( z ) = Bz 


( 11 ) 


Finally, the first term of Eq. 7 is written as 

J n [er'>(A)*, t n,i] dSt=-K-'f B t DB[z dU 


(12) 


The second term of Eq. 7 can also be derived in the same way to obtain a similar 
expression as in Eq. 12 as 



Bj'DBz dVL 


(13) 


Sip — 


Substituting Eqs. 12 and 13 into Eq. 7, the following expressions are obtained 

dz 


db t 


= j u - er*>(e)e^(A)V Jfc>Jfc ] dU 

JVe r 

= -K~ l y / B T DB'.z + Bj'DBz + B T DBzV k k dfl 

i 

= -K~ } y B t DB[ + Bj'DB + B T DBV k k dQ^j z 

= ~ k ~ a /;; /;; / / \b t db[ + bj'db + b t dbv^] i j i 


(14) 


where | J | is the determinant of the Jacobian matrix J which transforms the undeformed 
configuration into the natural coordinate system. The constraint functional change is 
equal to dz/db , , since all displacement sensitivities are calculated and the design change 
6b, is chosen as a unit number. 

For the implicit differentiation approach of Eq. 2, the derivative of the global stiffness 
matrix can be evaluated at the elemental level, i.e., 


dz 

db. 


= -K~ 1 



= -K ~ 1 



( 15 ) 
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where Kf is calculated numerically in natural coordinates as 


Kf = J*' j ^ J+* [ Bp DB | J | +BfDB' | J | +B t DB \J\' t ] d£dr]d^ (16) 


Comparing Eqs. 14 and 15, one sees that both are equivalent if the following expression 
is valid. 

\J\'i=\J\V k ,k (17) 

To prove Eq. 17 is true, one should notice that the determinant of the Jacobian can 
be separated into two parts. The first contribution is from the deformed to undeformed 

configuration, denoted by | J |<j, which depends on design. The other is the contribution 

of transformation from the undeformed global to local or natural configurations, denoted 
by | J |, which is independent of design. The relationship is expressed in the following 
form 

| J |r=| J\d\J\ (18) 

where r denotes the deformed configuration. Differentiating Eq. 18 with respect to design, 
one obtains 

I J l'r=l J li | J | (IS) 

It was shown in Ref. 4 that | J |Jj= at r = 0, if the design change is assumed to be 

equal to a unit vector. Thus, Eq. 19 is identical to the form of Eq. 17, and the equivalence 
of Eqs. 14 and 15 is proved. 

Another way to prove the validity of Eq. 17 is to carry out the differentiation directly 
by the definition of the Jacobian. Consider a two-dimensional case as an example, the 
right side of Eq. 17 is obained as 


tit,- _ ( dx dy dxdy\ f dV x dV y \ 

J 1 ** “ \d£ dn dr,dt)\dx dy) 

dy dV x dx dV y dy dV x dx dV y 

= &n~dl + d(W ~ dl~drf ~ d^~d( 


where the velocitv V T and V v are defined as 


(20) 


V, = 


dx 

dbi 

v - ihi 

« ~ dbi 


( 21 ) 


Substituting Eq. 21 into Eq. 20, the following expression is obtained 


_ d (dxdy Bydx \ 

/>c ’ k ~ db t dy d£ dr]) 

=i j i; 


( 22 ) 


This simple calculation also verifies that the relationship of Eq. 17 is valid. Note that the 
design change <56, is assumed to be unity in Eq. 17. 
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Numerical Verification and Comparison 

In this section, the equivalence of the domain method and the implicit differentiation 
is verified through a simple example. The accuracy of design sensitivities will be examined 
and compared through the refinement of the finite element mesh. 

A simple two-dimensional thin plate is considered as an example. The finite element 
configuration, dimensions, material properties, loading condition, and design variable are 
shown in Fig. 1. Design variable b is chosen to move the upper traction free boundary. 
The load of 100 lb is applied parabolically at the right of the plate. 

An 8-noded two-dimensional plane stress isoparametric element is employed for analy- 
sis. The boundary stresses and strains that appear in Eq. 3 are computed by extrapolat- 
ing linearly from the stresses at Gauss points, where the optimal or the best approximate 
stresses are located. Numerical results for design sensitivity of point A in the Y-direction 
for lxl, 2x2, 3x3. 4x4, 5x5 and 6x6 meshes are shown in Table 1. 

In Table 1, column 1 represents different finite element meshes and column 2 the 
displacement of point A in the Y-direction for the initial design, b. Columns 3 and 4 
represent the displacement sensitivities at point A of Fig. 1, using the boundary method 
(BM) of Eq. 3 and the domain integral (DM) of Eq. 7, respectively, for different meshes. 
Column 5 has the results using the implicit differentiation approach (IDA) of Eq. 2. The 
derivative of the global stiffness matrix is carried out by differentiating the element stiffness 
matrix, analytically. 

Fig. 2 shows the same results as in Table 1. From Fig. 2 and Table 1, one observes that 
the displacements and the sensitivities for the implicit approach (IDA) do not change much 
after 3x3 finite element mesh. However, the design sensitivity for the boundary method 
of the variational approach (BM) is still increasing at the limit of mesh refinement. This 
implies that the boundary method (BM) is more sensitive to the finite element results, 
although it provides the analytical formulation for sensitivities. And one concludes that 
the boundary method of the variational approach tends to yield better gradient estimates, 
once a more acciirate analysis is used and better boundary stresses are obtained. The 
same conclusion is also found in Refs. 5 and 8. 

Comparing column 4 with 5, one sees that the domain method results (DM) are very 
close to those obtained from the analytical implicit differentiation approach (IDA). Clearer 
interpretation can be observed from Fig. 2, which plots the displacement and displace- 
ment sensitivity versus finite element mesh size. This numerical agreement verifies the 
equivalence of the two approaches. 

In Ref. 8. the boundary the method was shown to be theoretically equivalent to the 

implicit approach, however, they yield slightly different results numerically as also shown 
in Table 1 and Fig. 2. The difference results from different numerical schemes for these two 
approaches, i.e., one uses the boundary, and the other the domain information. If consistent 
numerical schemes are used for the domain method and the implicit approach as in this 
report, they are shown to be equivalent, not only theoretically but also numerically. 

The disadvantage of the domain method is in computational aspects. Numerical eval- 
uation of Eq. 7 is more complicated than evaluation of Eq. 3, because Eq. 7 requires 
integration over the entire domain, whereas Eq. 3 requires integration only over the vari- 
able boundary. In addition, a velocity field which satisfies regularity properties must be 
defined in the domain. The choice of velocity for the physical domain is more difficult than 
that for the variable boundary. Although, a boundary layer scheme 9 and a displacement- 
like velocity field 10 were proposed to alleviate these problems, the domain method still 
requires more analyst and computational efforts. 
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Summary 

It is shown that accurate finite element analysis results in accurate design sensitivities. 
For the boundary method of the material derivative approach to shape design sensitivities, 
the accuracy of the finite element is more crucial, since the finite element method usually 
does not give accurate stresses on the boundary. 

The domain method is generally more accurate than the boundary method in the 
material derivative approach for evaluating the design sensitivities; however, a velocity 
field for the physical domain needs to be defined. The necessity of defining a domain 
velocity field and integrating the domain integral instead of the boundary integral, as in 
boundary method, requires both more analyst time and computational time. 

It is also shown that the domain method is equivalent, under certain assumptions, 
to the implicit differentiation approach not only theoretically but also numerically. The 
numerical equivalence is valid only if the numerical methods used for both approaches are 
consistent. 
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Table 1. Comparison of Design Sensitivity Accuracy 




BM 

DM 

IDA 

mesh 

displacement 

dv/db 

dv/db 

dv/db 

lxl 

2.495E-5 

-4.196E-6 

-4.843E-6 

-5.248E-6 

2x2 

2.760E-5 

-4.518E-6 

-5.167E-6 

-5.235E-6 

3x3 

2.841E-5 

-4.856E-6 

-5.369E-6 

-5.375E-6 

4x4 

2.845E-5 

-4.995E-6 

-5.391E-6 

-5.394E-6 

5x5 

2.854E-5 

-5.093E-6 

-5.412E-6 

-5.413E-6 

6x6 

2.856E-5 

-5.158E-6 

-5.425E-6 

-5.426E-6 




finite element meshes 
(8-node plane stress element) 


Y 


Fig. 1 Square Thin Plate 







